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We use quantum Monte Carlo simulations to determine the finite temperature phase diagram and 
to investigate the thermal and quantum melting of stripe phases in a two-dimensional hard-core 
boson model. At half filling and low temperatures the stripes melt at a first order transition. In the 
doped system, the melting transitions of the smectic phase at high temperatures and the superfiuid 
smectic (supersolid) phase at low temperatures are either very weakly first order, or of second order 
with no clear indications for an intermediate nematic phase. 
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Stripe phases of lattice models with broken rotational 
and translational symmetry can melt in two qualitatively 
different ways. One scenario is that both symmetries can 
be restored at a single first or second order transition, and 
the stripes melt directly into a normal fluid or superfluid 
phase. The other scenario is that first the translational 
symmetry is restored when the striped solid melts into 
a nematic (liquid crystal) phase with broken rotational 
symmetry. The rotational symmetry is then restored in 
a second melting transition of the nematic phase. 

This quantum lattice problem shows similarities to 
the long standing problem of the melting of a two di- 
mensional crystal into a continuum model, where there 
is also either a first order melting or the Kosterlitz- 
Thouless-Halperin-Nelson- Young (KTHNY) scenario of 
two Kosterlitz-Thouless transitions with an intervening 
hexactic phase Ij. Clear results for the classical version 
of this continuum problem were obtained only recently 
in a simple model of hard disks . 

Current interest in quantum mechanical stripe phases 
and their melting stems from the experimental observa- 
tion of stripes in some high- Tc superconductors Q and 
from the question whether and how they are related to 
the occurrence of high temperature superconductivity. 
Numerically, stripe phases have been found to be com- 
petitive ground states of t-J-\ike models 130. Analyti- 
cally some theories of high temperature superconductiv- 
ity are closely linked to the existence of stripe and ne- 
matic phases 0. While it is hard to study stripe phases 
in strongly correlated fermionic models, because of the 
negative sign problem of quantum Monte Carlo, we can 
more accurately investigate bosonic stripes using mod- 
ern quantum Monte Carlo algorithms • Such bosonic 
models can appear as effective low energy models neglect- 
ing nodal quasiparticles 0, 0|. Like the cuprates, these 
bosonic models show competition and in some models 
coexistence of superfluidity and charge order. 

In this Letter we focus on the simplest bosonic quan- 
tum model exhibiting stripe order and determine its finite 
temperature phase diagram. We carefully investigate its 
thermal and quantum melting transitions to address the 
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FIG. 1: Ground-state phase diagram of the hard-core boson 
Hubbard model Eq. Q as function of t/V2 and /x/Vg. Due 
to particle-hole symmetry, the phase diagram is symmetric 
around the half filling line (density p — 1/2 at /i = 2V2) and 
the lower half is shown. The thick lines 1) and 2) indicate 
the cuts along which we show the finite-temperature phase 
diagrams in Fig. |21 

question of how stripe melting occurs in a quantum model 
of bosonic stripes. The particular hard-core boson Hub- 
bard Hamiltonian we study on a two-dimensional square 
lattice is 

H — ^(ajaj + ajfli) + V2 ^ nirij — /i^nj, (1) 

(iJ> {(iJ» i 

where O; (oi) is the creation (annihilation) operator for 
hard-core bosons, ni = a[ai the number operator, V2 > 
the next nearest neighbor Coulomb repulsions and /i the 
chemical potential. The Hamiltonian ^ is equivalent to 
an anisotropic spin-1/2 model with a nearest neighbor 
Xy-like interaction J^y — 2t and next nearest neighbor 
Ising interaction = V2 in a magnetic field h = 2V2 — /i. 

In Fig. n we review the ground-state phase diagram of 
this model which was previously studied by mean-field, 
renormalization g roup and local-update quantum Monte 
Carlo simulations [iCtlllL 12] . At half filling (density p — 
1/2 at /i = 2V2) the ground state is a smectic for t/V2 < 
0.446 ± 0.006 and the quantum melting transition at low 
temperatures was found to be of first order (translation 
symmetry is broken only in one dimension, perpendicular 
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FIG. 3: Histograms of the stripe order parameter Os clearly 
show the double-peak structure of a first order transition 
which gets more pronounced as the system size is increased. 



FIG. 2: (a) Finite temperature phase diagram as function of 
t/V2 of the half filled model and (b) as function of the density 
p along the respective lines of Fig. The normal fluid and 
superfluid phases are denoted by the symbols NF and SF. 



to the stripes, hence this phase is a smectic) . 

Doping this smectic stripe crystal a stable "supersolid" 
phase with coexisting stripe order and superfluidity was 
found, in which the vacancies doped into the stripes form 
a superfluid. 

Before going into details we summarize our key re- 
sults by presenting the finite temperature phase diagrams 
along the lines indicated in the ground state phase dia- 
gram Fig. n At half filling (see Fig. we find that 
the transition between smectic and superfluid is of first 
order at all temperatures. At higher temperatures the 
transition between smectic and normal fluid is of first or- 
der around t/V2 ~ 0.5, but then becomes very weakly 
first order, before turning second order as t/V2 0. In 
the doped system (see Fig. I^d) coexistence between the 
smectic phase and superfiuidity can be observed at low 
temperatures around half filling. Here the phase transi- 
tions are all either very weakly first order or of second 
order. Rotational and translational symmetry breaking 
occur, within the accuracy of our simulations, at the same 
point. 

To determine the finite-temperature phase diagram we 
have used quantum Monte Carlo simulations, employing 
a directed loop quantum Monte Carlo algorithm in the in 
the stochastic series expansion (SSE) representation 0- 
The simulations were performed in the grand canonical 
ensemble and do not suffer from any systematic errors 
apart from finite size effects. In contrast to the loop 
algorithm 1^ the directed loop algorithm remains effec- 
tive away from half filling where the loop algorithm slows 
down exponentially. 

We determine stripe order by measuring the order pa- 
rameter 



Os 



5„(^,0)-H5n(0,^) 



(2) 



wave vector {kx,ky). To investigate a nematic phase, 
we have to look for rotational symmetry breaking in the 
kinetic energy or in the local charge correlations, using 
as order parameters 



Oh — — at ^a/ , , N — at ^a, ,■,■,+ H.c. (3) 

^ y (x,y) (x + l.,y) (x,y) (x.,y+l) V 

{x,y) 

or alternatively 



On = ^ «(a;,y)»^(K+l,a) " '^(£c,y)'^(a;,y-|-l) (4) 
{x,y) 



where n 



{x,y) 



0^{x,y)0'{x,y) 



is the boson number opera- 



tor at lattice site {x^y) and H.c. denotes the Hermitian 
conjugate. 

The extent of the superfiuid phase is determined by 
measuring the superfiuid (number) density 



Ps = mT{W^) 



(5) 



where Sn{kx,ky) is the charge structure factor at the 



where W is the winding number in one of the directions 
and m the mass of a boson, ps is finite in the superfiuid 
phase and jumps to zero from a universal value -^rnTc at 
the Kosterlitz-Thouless transition. 

We now discuss the phase diagrams and the nature 
of the phase transitions in more detail, starting with the 
half filled system at ^ = 2V2, where a first order quantum 
phase transition was found at ¥2/1 — 2.24±0.03, improv- 
ing the previous estimate • When raising the temper- 
ature we find that the transition remains of first order as 
the histograms for both the stripe order parameter Os 
and the rotational symmetry breaking order parameter 
On show two clearly separated peaks corresponding to 
the two coexisting phases (see Fig. 13 . Careful finite 
size scaling of the peak distances on system sizes L = 8 
up to L = 32 show that the two peaks stay separate as 
expected for a first order transition and do not merge 
as in a second order one. For both order parameters we 
find the same phase transition line, confirming a melt- 
ing in a single strong first order transition up to at least 
T/t = 2/3 [from point (D) to (B) in Fig.[2li)]. 
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FIG. 4: Fourth order cumulants for Os and On as function 
of the temperature (Vi = 0, V2 = 3t). The dashed horizontal 
line shows the critical value for a second order transition in the 
Ising universality class. The crossing point of the 4-th order 
cumulant ratios for Os is not consistent with this value. 



To investigate whether the stripes melt into a super- 
fluid or normal fluid we measure the superfluid density 
Ps at the coexistence line for configurations in the fluid 
phase (determined by their value of Os)- At a tempera- 
ture T/t = 2/3 [point (B)], ps drops below the universal 
jump value ^rnT as a function of system size L and will 
thus scale to zero in the infinite system. Hence the transi- 
tion at this temperature is of first order to a normal fluid 
phase in contrast to T = 0.5t [point (C)] where ps is well 
above the universal value for all system sizes indicating 
a direct first order transition between the smcctic and 
a superfluid phase |l3| . There is thus a tricritical point 
between points B and C where the Kosterlitz-Thouless 
phase transition of the superfluid becomes first order. 

At higher temperatures and larger repulsion V2 it be- 
comes harder to determine the order of the phase transi- 
tion. At V2 = 3t no double-peak structure could be seen 
in the histograms for lattice sizes up to i = 32 and the 
transition is thus either second order or a very weak first 
order transition. In the limit of infinite repulsion V2 the 
lattice decouples into two sub-lattices, each of which is 
equivalent to an Ising antiferromagnet with a second or- 
der melting transition. We thus explore the possibility of 
a second order transition. For a second order transition 
the 4-th order cumulant ratios C4 = 1— (C*) /3(0^)^ have 
a size independent value at the transition point, which 
can be found as the crossing points of the cumulants for 
different system sizes. Our results (see Fig. 0J show that 
we indeed observe such crossings, which is and indication 
for second order transition. Both the 4-th order cumulant 
ratios for both On and Os cross at the same tempera- 
ture within the accuracy of our results, indicating that 
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FIG. 5: Superfluid density parallel {Ps,\{) and perpendicular 
{ps,±) to the stripes for L — 8, V^jt — 5 and T = t/6. 



translational symmetry breaking (measured by Os) and 
rotational symmetry breaking (measured by O^r) happen 
at the same or at very close temperatures. 

The phase diagram of the doped system away from half 
filling shows an additional "supersolid" phase where the 
vacancies doped into the stripes form a superfluid smec- 
tic with broken translational and rotational symmetry 
as well as a finite superfluid density. Since the rota- 
tional symmetry is spontaneously broken in the smectic 
phase, the superfluid becomes anisotropic and Eq. © 
for the superfluid density needs to be modified. The 
geometric mean of the superfluid densities perpendicu- 
lar and parallel to the stripes ps = \/Ps,\\ ' Ps,i^ is the 
quantity exhibiting the universal jump at the Kosterlitz- 
Thouless transition and was used in the determination of 
the phase transition. Care was taken that the anisotropy 
did not become too large to introduce systematic errors 
as pointed out by Prokof'ev and Svistunov 0]. The su- 
perfluid densities parallel {Ps,\\) and perpendicular (ps._l) 
to the stripe order were measured for each configuration 
from the winding numbers in either x or y direction de- 
pending on whether S'„(0, tt) was larger or smaller than 
^„(^,0). 

Although superfluidity coexists with smectic order, it 
is strongly suppressed and Tc and p^ vanish linearly 
as half filling is approached. Using = l/2t and 
wj_ = 2t'i /(4V2) ^^'^ boson masses parallel and per- 
pendicular to the stripes, respectively, we obtain p^ — 
(0.78 ± 0.07) |p — i| consistent with the value obtained 
previously for a dilute gas of bosons . 

In order to investigate the nature of the doping-driven 
melting of the smectic phase we repeat a similar pro- 
cedure as at half filling. In contrast to the half filled 
case, with a clear jump in the values of Os and On at 
the phase transition, or in the doped case of the nearest 
neighbor model with a jurnp in density p at the doping- 
driven melting transition [l6j, no such jump or double- 
peak structure was observed here. We checked the order 
parameters Os, On, Ok, and the density p in systems 
with size up to i = 48. Repeating the simulations of 
Ref. where a local update algorithm was used in- 
stead of our nonlocal algorithm we cannot reproduce the 



4 



jump in the superfluid density seen there. Instead we 
find the smooth behavior shown in Fig. ]E\ for the same 
system, which is also in agreement with analytical consid- 
erations using mean-field and spin-wave analysis or 
renormalization group calculations ^3 ■ Since our results 
are well converged and do not suffer from any systematic 
error, we believe the differences to be due to the method 
used to calculate the superfluid densities in Ref. 0|. 

At T = t/6 we checked whether rotational and trans- 
lational symmetry breaking occur at the same critical 
density by comparing estimates obtained using Binder 
cumulant ratios for Os and On as well as a determina- 
tion of the critical point from the maxima of the suscep- 
tibility associated with Ok- The estimates agree within 
the error bars: pc = 0.2468 ± 0.0025 at T = t/6. 

Finally we want to discuss our results in view of the 
three possibilities for the nature of the melting transi- 
tions of i) a first order transition, ii) a single second 
order transition or iii) two separate second order tran- 
sitions. Close to half filling and at low temperatures we 
have very clear evidence for a single first order transition. 
The results for thermal melting at higher temperatures 
and for the doping-driven melting at low temperatures 
are not as clear-cut. While a very weak first order tran- 
sition is always possible, it is more likely that we have 
second order transitions, especially in view of the fact 
that in the limit t/V2 — > oo there will be a second order 
transition in the Ising universality class. One way to dis- 
tinguish between the scenarios (ii) of one single transition 
or (iii) of two second order transitions with an intervening 
nematic phase is that in the latter case both transitions 
should be in the Ising universality class since they each 
break a Z2 symmetry. At the transition we would thus 
expect that the Binder cumulant ratios at Tc takes on the 
value C4 = 0.6106900(1) The fact that our results 

for the cumulant ratio of Os are close but not identical 
to this universal value for an Ising transition rather in- 
dicates secnario (ii), a single second order transition in 
a different universality class. In the doped system our 
data both for the cumulant ratios and the scaling of the 
susceptibility of Ok are actually closer to the three dimen- 
sioonal XY universality class predicted by Ref. than 
the Ising universality class. More extensive simulations 
on larger systems, using newly developed flat histogram 
methods for quantum systems |l8| will be needed to more 
accurately determine the universality class of this phase 
transition. 

While the exact order and universality class of the 
melting transition are hard to determine - which is not 
surprising given the difficulties known from the similar 
classical problem in the continuum ~ our results are in- 
teresting also with respect to the existence of a nematic 
phase. Such a phase must be restricted to a very narrow 
temperature and doping regime, smaller than the resolu- 
tion of our simulations. It might be possible to stabilize a 



nematic phase with additional terms in the Hamiltonoan. 
One suggestion 0| is to add a term proportional to 
the square of the nermatic order parameter —VO^, to 
stabilize a nematic phase. This term gives two con- 
tributions: a nearest neighbor repulsion 2y^^jj^ninj 
and an additional next nearest neighbor hopping term 
^J2{{i j}) (^i + aj ai)- These frustrating terms, the sec- 
ond of which unfortunately causes a negative sign prob- 
lem for quantum Monte Carlo simulations, might be im- 
portant for the stability of an extended nematic phase. 
It is interesting to compare this model with the related 
frustrated Heisenberg model on a square lattice which 
however exhibits translational and rotational symmetry 
breaking only in the ground state at T = |2(il |. 
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